function varargout=compute_MBC_shock10(varargin)
% This function computes MBC shock (Main Business Cycle) denoted as q
% input
% Compared to compute_MBC_shock, this allows multiple variables 
% varargin{1}=out: VAR result
% varargin{2}=ps: parameters
% varargin{3}=opt: option
% 
% output is rs including q
% 
% 2021/10/2 copied from compute_MBC_shock6
% 
%

%--------------------------------------------------------------------------
% input
out=varargin{1};
ps=varargin{2};
opt=varargin{3};
%rs=varargin{4};
%--------------------------------------------------------------------------
% compute H 
H=zeros(ps.N,ps.N); % integral of spectral density D(w) associated with q
D_vec=zeros(ps.ngrid,ps.N,ps.N); % D(w)
D_all_vec=zeros(ps.ngrid,1); % spectral density at each w

qDq_vec=zeros(ps.ngrid,ps.N,ps.N);

bin_vec=[];%zeros(ps.nbin,1); % bin values, not used anymore


for i=1:ps.ngrid % grids [0,2*pi]

    
% compute D(w)    
% tmp=compute_Dw6(ps.grid(i),out,ps,opt,rs);
tmp=compute_Dw10(ps.grid(i),out,ps,opt);

% to compute q'D(w)q (spectral density associated with q at w)
qDq_vec(i,:,:)=tmp.D; % D(w) at w

% q shock, matrix D(w) in [w1,w2] with width
D_vec(i,:,:)=tmp.D*ps.igrid(i)*(2*pi/ps.ngrid); % D(w)*width,0 if grid outside [w1,w2]

% all shocks (scalar, spectral density at w)
D_all_vec(i)=tmp.D_all;

H=H+squeeze(D_vec(i,:,:)); % integral of D(w) in [w1,w2]

end


% compute eigenvector

% only largest eigenvalue, faster
% [eig_vec,eig_val]=eigs(H,1); % 1st largest eigenvalue and eigenvector q 
% q=eig_vec(:,1); % 1st largest eigenvector

% if want to get all eigenvalues, slower
[eig_vec,eig_val]=eig(H); % all eigenvalues
[d,ind]=sort(diag((eig_val)),'descend'); % largest to smallest order
eig_val=eig_val(ind,ind); % reorder eigenvalue
eig_vec=eig_vec(:,ind); % reorder eigenvector
q=eig_vec(:,1); % 1st largest eigenvalue 
%eig_val=[];
%eig_vec=[];

% normalize the eigenvector
q=q/norm(q); % |q|=1 

% if q(ps.k)<0 % q(k)>0 : k-th row of q is positive
%     q=-q;
% end
S_tilda=chol(out.SIGMA,'lower');
tmp=S_tilda*q;

if tmp(ps.k2(1))<0
    q=-q;
end

q1=q;

% compute D(w) for q
qDq=zeros(ps.ngrid,1);
for i=1:ps.ngrid
    qDq(i)=q'*squeeze(qDq_vec(i,:,:))*q;
end

% Objective function
%GAMMASTAR_CURR = @(x) -GAMMASTAR(H,B,A_0,SIGMA,x);
obj_fun = @(x) -obj_fun_ex_q2(x,H,ps,opt);%,rs);

%ub = [0;repmat(Inf,[k-1 1])];
%lb = [0;repmat(-Inf,[k-1 1])];
A_eq=opt.q1';
%A_eq=A_eq(2:end); % eliminate the first element
b_eq=0;
%options = optimset('Algorithm','active-set','Display','iter');
options = optimset('Algorithm','active-set','Display','off');

[x,fval,exitflag,output] = fmincon(obj_fun,...
    sqrt(1/(ps.N))*ones(ps.N,1),... % inital value
    [],[],... % inequality constraints
    A_eq,b_eq,... % equality constraints
    [],[],... % lower and uppear abounds
    @(x) mycon(x),... % nonlinear constraints [c<=0, ceq=0] 
    options);

q2 = x;


% normalize the eigenvector
q2=q2/norm(q2); % |q|=1 

% if q(ps.k)<0 % q(k)>0 : k-th row of q is positive
%     q=-q;
% end
S_tilda=chol(out.SIGMA,'lower');
tmp=S_tilda*q2;

if tmp(ps.k2(1))<0
    q2=-q2;
end


%--------------------------------------------------------------------------
% output
rs.q1=q1;
rs.q2=q2;
rs.eig_vec=eig_vec;
rs.eig_val=eig_val;
rs.D_all_vec=D_all_vec;
rs.bin_vec=bin_vec;
rs.qDq=qDq;
varargout{1}=rs;




